/*-
 * Copyright (c) 2013 Shang-Rong Cai. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHORS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 * $ MTBLAS: stblas_2_gmv.c    Rev. 0.0.25    2013-10-23T17:16+08:00 $
 */

#include <stdlib.h>
#ifdef TRACE_INFO
#include <stdio.h>
#endif

#include <bscerr.h>

#include <exec_isa.h>
#include <cor_ctrl.h>
#include <stblas.h>

/* Locally Common Function For Setting Operands of Matrix-Vector Multiplication */
void stgmv_opchunks_mv(int flag, step_t *step_one, step_t *step_fix, step_t *step_lda, step_t *step_incx, addr_t ma, addr_t vx, addr_t buf, size_t c_sz, opchunk_t *pk_a, opchunk_t *pk_b, opchunk_t *pk_c)
{
    pk_a->sz = c_sz;
    pk_a->ra.addr = ma;

    pk_a->rb.addr = vx;
    pk_a->rb.step = *step_incx;
    pk_a->rc.addr = buf;
    pk_a->rc.step = *step_fix;
    pk_a->rd.addr = NULL;
    pk_a->rd.step = *step_fix;

    pk_b->sz = c_sz;
    pk_b->ra.addr = ma;

    pk_b->rb.addr = vx;
    pk_b->rb.step = *step_fix;
    pk_b->rc.addr = buf;
    pk_b->rc.step = *step_one;
    pk_b->rd.addr = NULL;
    pk_b->rd.step = *step_fix;

    pk_c->ra.step = *step_fix;
    pk_c->rb.addr = step_incx;
    pk_c->rb.step = *step_fix;
    pk_c->rc.addr = step_fix;
    pk_c->rc.step = *step_fix;
    pk_c->rd.addr = step_fix;
    pk_c->rd.step = *step_fix;

    if (flag == GEMV_CATEGORY_0)
    {
        pk_a->ra.step = *step_one;
        pk_b->ra.step = *step_lda;
        pk_c->ra.addr = step_one;
    }
    else
    {
        pk_a->ra.step = *step_lda;
        pk_b->ra.step = *step_one;
        pk_c->ra.addr = step_lda;
    }
}


int stblas_gemvf(const int rdr, const int trn, const size_t rz, const size_t cz, const float alpha, const float *ma, const size_t lda, const float *vx, const step_t incx, const float beta, float *vy, const step_t incy)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(float);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(float) * lda;
    const step_t step_incx = sizeof(float) * incx;

    const float buf_init = 0;

    size_t r_sz;  /* Length of result vector */
    size_t c_sz;  /* Length of FMAC */

    int flag;  /* Type of matrix order and transpose */

    float *buf;  /* Buffer vector for storing ma * vx */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to modifier of address */
    operand_t rc;  /* Usual operand or point to modifier of step */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_gemvf()\n");
#endif

    /* Set Controls */
    ierr = cor_ctrl_gemv_flag(rdr, trn, &flag);

    if (ierr != ERRNO_UNDEFINED)
    {
        goto EXIT_STBLAS_GEMVF;
    }

    ierr = cor_ctrl_gemv_size(trn, rz, cz, &r_sz, &c_sz);

    /* Allocate Memory For Buffer: buf[] */
    buf = (float*) calloc(r_sz, sizeof(float));

    /* buf[] = 0 */
    ra.addr = (void*) &buf_init;
    ra.step = 0;
    rb.addr = (void*) buf;
    rb.step = sizeof(float);

    execute((void*) op2_copy_ff, r_sz, &ra, &rb, NULL, NULL);

    /* buf[] = ma * vx */
    pk_a.oc = (addr_t) op3_fmac_ff;
    pk_b.oc = (addr_t) op3_fmac_ff;

    stgmv_opchunks_mv(flag, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (step_t*) &step_incx, (addr_t) ma, (addr_t) vx, (addr_t) buf, c_sz, &pk_a, &pk_b, &pk_c);

    ra.addr = (addr_t) &pk_a;
    ra.step = step_fix;
    rb.addr = (addr_t) &pk_b;
    rb.step = step_fix;
    rc.addr = (addr_t) &pk_c;
    rc.step = step_fix;

    execute((void*) op4_exec_op, r_sz, &ra, &ra, &rb, &rc);

    /* beta * vy */
    ra.addr = (void*) &beta;
    ra.step = 0;
    rb.addr = (void*) vy;
    rb.step = sizeof(float) * incy;
    rc.addr = (void*) vy;
    rc.step = sizeof(float) * incy;

    execute((void*) op3_mul_ff, r_sz, &ra, &rb, &rc, NULL);

    /* vy += alpha * buf[] */
    ra.addr = (void*) &alpha;
    ra.step = 0;
    rb.addr = (void*) buf;
    rb.step = sizeof(float);
    rc.addr = (void*) vy;
    rc.step = sizeof(float) * incy;
    execute((void*) op3_fmac_ff, r_sz, &ra, &rb, &rc, NULL);

    free(buf);

    EXIT_STBLAS_GEMVF:
#ifdef TRACE_INFO
    printf("    stblas_gemvf() Passed!\n");
#endif

    return ierr;
}


int stblas_gemv(const int rdr, const int trn, const size_t rz, const size_t cz, const double alpha, const double *ma, const size_t lda, const double *vx, const step_t incx, const double beta, double *vy, const step_t incy)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(double);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(double) * lda;
    const step_t step_incx = sizeof(double) * incx;

    const double buf_init = 0;

    size_t r_sz;  /* Length of result vector */
    size_t c_sz;  /* Length of FMAC */

    int flag;  /* Type of matrix order and transpose */

    double *buf;  /* Buffer vector for storing ma * vx */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to modifier of address */
    operand_t rc;  /* Usual operand or point to modifier of step */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_gemv()\n");
#endif

    /* Set Controls */
    ierr = cor_ctrl_gemv_flag(rdr, trn, &flag);

    if (ierr != ERRNO_UNDEFINED)
    {
        goto EXIT_STBLAS_GEMV;
    }

    ierr = cor_ctrl_gemv_size(trn, rz, cz, &r_sz, &c_sz);

    /* Allocate Memory For Buffer: buf[] */
    buf = (double*) calloc(r_sz, sizeof(double));

    /* buf[] = 0 */
    ra.addr = (void*) &buf_init;
    ra.step = 0;
    rb.addr = (void*) buf;
    rb.step = sizeof(double);

    execute((void*) op2_copy_df, r_sz, &ra, &rb, NULL, NULL);

    /* buf[] = ma * vx */
    pk_a.oc = (addr_t) op3_fmac_df;
    pk_b.oc = (addr_t) op3_fmac_df;

    stgmv_opchunks_mv(flag, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (step_t*) &step_incx, (addr_t) ma, (addr_t) vx, (addr_t) buf, c_sz, &pk_a, &pk_b, &pk_c);

    ra.addr = (addr_t) &pk_a;
    ra.step = step_fix;
    rb.addr = (addr_t) &pk_b;
    rb.step = step_fix;
    rc.addr = (addr_t) &pk_c;
    rc.step = step_fix;

    execute((void*) op4_exec_op, r_sz, &ra, &ra, &rb, &rc);

    /* beta * vy */
    ra.addr = (void*) &beta;
    ra.step = 0;
    rb.addr = (void*) vy;
    rb.step = sizeof(double) * incy;
    rc.addr = (void*) vy;
    rc.step = sizeof(double) * incy;

    execute((void*) op3_mul_df, r_sz, &ra, &rb, &rc, NULL);

    /* vy += alpha * buf[] */
    ra.addr = (void*) &alpha;
    ra.step = 0;
    rb.addr = (void*) buf;
    rb.step = sizeof(double);
    rc.addr = (void*) vy;
    rc.step = sizeof(double) * incy;
    execute((void*) op3_fmac_df, r_sz, &ra, &rb, &rc, NULL);

    free(buf);

    EXIT_STBLAS_GEMV:
#ifdef TRACE_INFO
    printf("    stblas_gemv() Passed!\n");
#endif

    return ierr;
}


int stblas_gemvl(const int rdr, const int trn, const size_t rz, const size_t cz, const long double alpha, const long double *ma, const size_t lda, const long double *vx, const step_t incx, const long double beta, long double *vy, const step_t incy)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(long double);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(long double) * lda;
    const step_t step_incx = sizeof(long double) * incx;

    const long double buf_init = 0;

    size_t r_sz;  /* Length of result vector */
    size_t c_sz;  /* Length of FMAC */

    int flag;  /* Type of matrix order and transpose */

    long double *buf;  /* Buffer vector for storing ma * vx */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to modifier of address */
    operand_t rc;  /* Usual operand or point to modifier of step */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_gemvl()\n");
#endif

    /* Set Controls */
    ierr = cor_ctrl_gemv_flag(rdr, trn, &flag);

    if (ierr != ERRNO_UNDEFINED)
    {
        goto EXIT_STBLAS_GEMVL;
    }

    ierr = cor_ctrl_gemv_size(trn, rz, cz, &r_sz, &c_sz);

    /* Allocate Memory For Buffer: buf[] */
    buf = (long double*) calloc(r_sz, sizeof(long double));

    /* buf[] = 0 */
    ra.addr = (void*) &buf_init;
    ra.step = 0;
    rb.addr = (void*) buf;
    rb.step = sizeof(long double);

    execute((void*) op2_copy_lf, r_sz, &ra, &rb, NULL, NULL);

    /* buf[] = ma * vx */
    pk_a.oc = (addr_t) op3_fmac_lf;
    pk_b.oc = (addr_t) op3_fmac_lf;

    stgmv_opchunks_mv(flag, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (step_t*) &step_incx, (addr_t) ma, (addr_t) vx, (addr_t) buf, c_sz, &pk_a, &pk_b, &pk_c);

    ra.addr = (addr_t) &pk_a;
    ra.step = step_fix;
    rb.addr = (addr_t) &pk_b;
    rb.step = step_fix;
    rc.addr = (addr_t) &pk_c;
    rc.step = step_fix;

    execute((void*) op4_exec_op, r_sz, &ra, &ra, &rb, &rc);

    /* beta * vy */
    ra.addr = (void*) &beta;
    ra.step = 0;
    rb.addr = (void*) vy;
    rb.step = sizeof(long double) * incy;
    rc.addr = (void*) vy;
    rc.step = sizeof(long double) * incy;

    execute((void*) op3_mul_lf, r_sz, &ra, &rb, &rc, NULL);

    /* vy += alpha * buf[] */
    ra.addr = (void*) &alpha;
    ra.step = 0;
    rb.addr = (void*) buf;
    rb.step = sizeof(long double);
    rc.addr = (void*) vy;
    rc.step = sizeof(long double) * incy;
    execute((void*) op3_fmac_lf, r_sz, &ra, &rb, &rc, NULL);

    free(buf);

    EXIT_STBLAS_GEMVL:
#ifdef TRACE_INFO
    printf("    stblas_gemvl() Passed!\n");
#endif

    return ierr;
}


int stblas_cgemvf(const int rdr, const int trn, const size_t rz, const size_t cz, const float complex alpha, const float complex *ma, const size_t lda, const float complex *vx, const step_t incx, const float complex beta, float complex *vy, const step_t incy)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(float complex);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(float complex) * lda;
    const step_t step_incx = sizeof(float complex) * incx;

    const float complex buf_init = 0;

    size_t r_sz;  /* Length of result vector */
    size_t c_sz;  /* Length of FMAC */

    int flag;  /* Type of matrix order and transpose */

    float complex *buf;  /* Buffer vector for storing ma * vx */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to modifier of address */
    operand_t rc;  /* Usual operand or point to modifier of step */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_cgemvf()\n");
#endif

    /* Set Controls */
    ierr = cor_ctrl_gemv_flag(rdr, trn, &flag);

    if (ierr != ERRNO_UNDEFINED)
    {
        goto EXIT_STBLAS_CGEMVF;
    }

    ierr = cor_ctrl_gemv_size(trn, rz, cz, &r_sz, &c_sz);

    /* Allocate Memory For Buffer: buf[] */
    buf = (float complex*) calloc(r_sz, sizeof(float complex));

    /* buf[] = 0 */
    ra.addr = (void*) &buf_init;
    ra.step = 0;
    rb.addr = (void*) buf;
    rb.step = sizeof(float complex);

    execute((void*) op2_copy_fc, r_sz, &ra, &rb, NULL, NULL);

    /* buf[] = ma * vx */
    if (trn == MTBLAS_TRANS_CONJ)
    {
        pk_a.oc = (addr_t) op3_conjfmac_fc;
        pk_b.oc = (addr_t) op3_conjfmac_fc;
    }
    else
    {
        pk_a.oc = (addr_t) op3_fmac_fc;
        pk_b.oc = (addr_t) op3_fmac_fc;
    }

    stgmv_opchunks_mv(flag, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (step_t*) &step_incx, (addr_t) ma, (addr_t) vx, (addr_t) buf, c_sz, &pk_a, &pk_b, &pk_c);

    ra.addr = (addr_t) &pk_a;
    ra.step = step_fix;
    rb.addr = (addr_t) &pk_b;
    rb.step = step_fix;
    rc.addr = (addr_t) &pk_c;
    rc.step = step_fix;

    execute((void*) op4_exec_op, r_sz, &ra, &ra, &rb, &rc);

    /* beta * vy */
    ra.addr = (void*) &beta;
    ra.step = 0;
    rb.addr = (void*) vy;
    rb.step = sizeof(float complex) * incy;
    rc.addr = (void*) vy;
    rc.step = sizeof(float complex) * incy;

    execute((void*) op3_mul_fc, r_sz, &ra, &rb, &rc, NULL);

    /* vy += alpha * buf[] */
    ra.addr = (void*) &alpha;
    ra.step = 0;
    rb.addr = (void*) buf;
    rb.step = sizeof(float complex);
    rc.addr = (void*) vy;
    rc.step = sizeof(float complex) * incy;
    execute((void*) op3_fmac_fc, r_sz, &ra, &rb, &rc, NULL);

    free(buf);

    EXIT_STBLAS_CGEMVF:
#ifdef TRACE_INFO
    printf("    stblas_cgemvf() Passed!\n");
#endif

    return ierr;
}


int stblas_cgemv(const int rdr, const int trn, const size_t rz, const size_t cz, const double complex alpha, const double complex *ma, const size_t lda, const double complex *vx, const step_t incx, const double complex beta, double complex *vy, const step_t incy)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(double complex);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(double complex) * lda;
    const step_t step_incx = sizeof(double complex) * incx;

    const double complex buf_init = 0;

    size_t r_sz;  /* Length of result vector */
    size_t c_sz;  /* Length of FMAC */

    int flag;  /* Type of matrix order and transpose */

    double complex *buf;  /* Buffer vector for storing ma * vx */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to modifier of address */
    operand_t rc;  /* Usual operand or point to modifier of step */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_cgemv()\n");
#endif

    /* Set Controls */
    ierr = cor_ctrl_gemv_flag(rdr, trn, &flag);

    if (ierr != ERRNO_UNDEFINED)
    {
        goto EXIT_STBLAS_CGEMV;
    }

    ierr = cor_ctrl_gemv_size(trn, rz, cz, &r_sz, &c_sz);

    /* Allocate Memory For Buffer: buf[] */
    buf = (double complex*) calloc(r_sz, sizeof(double complex));

    /* buf[] = 0 */
    ra.addr = (void*) &buf_init;
    ra.step = 0;
    rb.addr = (void*) buf;
    rb.step = sizeof(double complex);

    execute((void*) op2_copy_dc, r_sz, &ra, &rb, NULL, NULL);

    /* buf[] = ma * vx */
    if (trn == MTBLAS_TRANS_CONJ)
    {
        pk_a.oc = (addr_t) op3_conjfmac_dc;
        pk_b.oc = (addr_t) op3_conjfmac_dc;
    }
    else
    {
        pk_a.oc = (addr_t) op3_fmac_dc;
        pk_b.oc = (addr_t) op3_fmac_dc;
    }

    stgmv_opchunks_mv(flag, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (step_t*) &step_incx, (addr_t) ma, (addr_t) vx, (addr_t) buf, c_sz, &pk_a, &pk_b, &pk_c);

    ra.addr = (addr_t) &pk_a;
    ra.step = step_fix;
    rb.addr = (addr_t) &pk_b;
    rb.step = step_fix;
    rc.addr = (addr_t) &pk_c;
    rc.step = step_fix;

    execute((void*) op4_exec_op, r_sz, &ra, &ra, &rb, &rc);

    /* beta * vy */
    ra.addr = (void*) &beta;
    ra.step = 0;
    rb.addr = (void*) vy;
    rb.step = sizeof(double complex) * incy;
    rc.addr = (void*) vy;
    rc.step = sizeof(double complex) * incy;

    execute((void*) op3_mul_dc, r_sz, &ra, &rb, &rc, NULL);

    /* vy += alpha * buf[] */
    ra.addr = (void*) &alpha;
    ra.step = 0;
    rb.addr = (void*) buf;
    rb.step = sizeof(double complex);
    rc.addr = (void*) vy;
    rc.step = sizeof(double complex) * incy;
    execute((void*) op3_fmac_dc, r_sz, &ra, &rb, &rc, NULL);

    free(buf);

    EXIT_STBLAS_CGEMV:
#ifdef TRACE_INFO
    printf("    stblas_cgemv() Passed!\n");
#endif

    return ierr;
}


int stblas_cgemvl(const int rdr, const int trn, const size_t rz, const size_t cz, const long double complex alpha, const long double complex *ma, const size_t lda, const long double complex *vx, const step_t incx, const long double complex beta, long double complex *vy, const step_t incy)
{
    int ierr = ERRNO_UNDEFINED;

    const step_t step_one = sizeof(long double complex);
    const step_t step_fix = 0;
    const step_t step_lda = sizeof(long double complex) * lda;
    const step_t step_incx = sizeof(long double complex) * incx;

    const long double complex buf_init = 0;

    size_t r_sz;  /* Length of result vector */
    size_t c_sz;  /* Length of FMAC */

    int flag;  /* Type of matrix order and transpose */

    long double complex *buf;  /* Buffer vector for storing ma * vx */

    operand_t ra;  /* Usual operand or point to opchunk to target of executing */
    operand_t rb;  /* Usual operand or point to modifier of address */
    operand_t rc;  /* Usual operand or point to modifier of step */

    opchunk_t pk_a;  /* Target of executing and modifying operations execution pack */
    opchunk_t pk_b;  /* Modifier of address */
    opchunk_t pk_c;  /* Modifier of step */

#ifdef TRACE_INFO
    printf("[F] stblas_cgemvl()\n");
#endif

    /* Set Controls */
    ierr = cor_ctrl_gemv_flag(rdr, trn, &flag);

    if (ierr != ERRNO_UNDEFINED)
    {
        goto EXIT_STBLAS_CGEMVL;
    }

    ierr = cor_ctrl_gemv_size(trn, rz, cz, &r_sz, &c_sz);

    /* Allocate Memory For Buffer: buf[] */
    buf = (long double complex*) calloc(r_sz, sizeof(long double complex));

    /* buf[] = 0 */
    ra.addr = (void*) &buf_init;
    ra.step = 0;
    rb.addr = (void*) buf;
    rb.step = sizeof(long double complex);

    execute((void*) op2_copy_lc, r_sz, &ra, &rb, NULL, NULL);

    /* buf[] = ma * vx */
    if (trn == MTBLAS_TRANS_CONJ)
    {
        pk_a.oc = (addr_t) op3_conjfmac_lc;
        pk_b.oc = (addr_t) op3_conjfmac_lc;
    }
    else
    {
        pk_a.oc = (addr_t) op3_fmac_lc;
        pk_b.oc = (addr_t) op3_fmac_lc;
    }

    stgmv_opchunks_mv(flag, (step_t*) &step_one, (step_t*) &step_fix, (step_t*) &step_lda, (step_t*) &step_incx, (addr_t) ma, (addr_t) vx, (addr_t) buf, c_sz, &pk_a, &pk_b, &pk_c);

    ra.addr = (addr_t) &pk_a;
    ra.step = step_fix;
    rb.addr = (addr_t) &pk_b;
    rb.step = step_fix;
    rc.addr = (addr_t) &pk_c;
    rc.step = step_fix;

    execute((void*) op4_exec_op, r_sz, &ra, &ra, &rb, &rc);

    /* beta * vy */
    ra.addr = (void*) &beta;
    ra.step = 0;
    rb.addr = (void*) vy;
    rb.step = sizeof(long double complex) * incy;
    rc.addr = (void*) vy;
    rc.step = sizeof(long double complex) * incy;

    execute((void*) op3_mul_lc, r_sz, &ra, &rb, &rc, NULL);

    /* vy += alpha * buf[] */
    ra.addr = (void*) &alpha;
    ra.step = 0;
    rb.addr = (void*) buf;
    rb.step = sizeof(long double complex);
    rc.addr = (void*) vy;
    rc.step = sizeof(long double complex) * incy;
    execute((void*) op3_fmac_lc, r_sz, &ra, &rb, &rc, NULL);

    free(buf);

    EXIT_STBLAS_CGEMVL:
#ifdef TRACE_INFO
    printf("    stblas_cgemvl() Passed!\n");
#endif

    return ierr;
}
